Gravitational Wavetrains in the Quasi-Equilibrium Approximation: 
A Model Problem in Scalar Gravitation 
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A quasi-equilibrium (QE) computational scheme was recently developed in general relativity to 
calculate the complete gravitational wavetrain emitted during the inspiral phase of compact binaries. 
The QE method exploits the fact that the the gravitational radiation inspiral timescale is much 
longer than the orbital period everywhere outside the ISCO. Here we demonstrate the validity and 
advantages of the QE scheme by solving a model problem in relativistic scalar gravitation theory. By 
adopting scalar gravitation, we are able to numerically track without approximation the damping of 
a simple, quasi-periodic radiating system (an oscillating spherical matter shell) to final equilibrium, 
and then use the exact numerical results to calibrate the QE approximation method. In particular, 
we calculate the emitted gravitational wavetrain three different ways: by integrating the exact 
coupled dynamical field and matter equations, by using the scalar-wave monopole approximation 
formula (corresponding to the quadrupole formula in general relativity), and by adopting the QE 
scheme. We find that the monopole formula works well for weak field cases, but fails when the fields 
become even moderately strong. By contrast, the QE scheme remains quite reliable for moderately 
strong fields, and begins to breakdown only for ultra-strong fields. The QE scheme thus provides 
a promising technique to construct the complete wavetrain from binary inspiral outside the ISCO, 
where the gravitational fields are strong, but where the computational resources required to follow 
the system for more than a few orbits by direct numerical integration of the exact equations are 
prohibitive. 



I. INTRODUCTION 

The construction of several new gravitational wave de- 
tectors, including the Laser Interferometer Gravitational 
Wave Observatory (LIGO), TAMA, VIRGO and GEO, 
may soon make gravitational wave astronomy a reality. 
The inspiral and coalescence of compact binaries, consist- 
ing of neutron stars or black holes, are among the most 
promising sources for detection by these observatories. It 
is expected that neutron star/neutron star binaries will 
spend approximately 16,000 cycles in the LIGO/ VIRGO 
frequency band, neutron star/black hole binaries about 
3,500, and black hole/black hole binaries about 600 [Q. 
To increase the likelihood of detection, and to extract in- 
formation from the signal, the binary inspiral has to be 
modeled theoretically and waveform templates have to 
be constructed. 

The evolution of compact binaries proceeds in different 
stages. By far the longest is the initial quasi-equilibrium 
inspiral stage, during which the compact objects move 
in nearly circular orbits, while the separation between 
them slowly decreases as energy is carried away by grav- 
itational radiation. The quasi-circular orbits become 
unstable at the innermost stable circular orbit (ISCO), 
where the inspiral enters a plunge and merger phase. 
The merger and coalescence happens on a dynamical 
timescale, and produces either a black hole or, for binary 
neutron stars, possibly a larger neutron star, which may 
collapse to a black hole at a later time. The final stage 
of the evolution is the ringdown phase, during which the 



merged object settles down to equilibrium. 

Two distinct approaches have commonly been em- 
ployed to analyze the adiabatic inspiral phase. Much 
progress has been made in post-Newtonian studies of 
compact binaries (see e.g. and references therein). 
Most of these approaches, however, approximate the 
compact objects as point sources, which neglects impor- 
tant finite-size effects which may be particularly impor- 
tant for neutron star binaries. Also, PN expansions may 
not converge sufficiently rapidly in the strong-field re- 
gion near the ISCO. Alternatively, compact binaries can 
be modeled numerically. Computational constraints cur- 
rently limit dynamical evolution calculations to at most 
a few orbits, so that there is no hope to simulate a com- 
plete inspiral. It is possible, however, to numerically 
model binaries in the adiabatic inspiral phase in a quasi- 
equilibrium (QE) approximation. 

The QE approximation is based on the assumption 
that the gravitational radiation reaction time scale is 
much longer than the orbital timescale, so that the bi- 
nary can be approximated to be in quasi-equilibrium (and 
quasi-circular orbit) on an orbital timescale. A similar 
approximation is routinely being used in stellar evolu- 
tion calculations. There, the evolutionary timescale is 
much longer than the hydrodynamical timescale, so that 
the star can safely be assumed to be in quasi-equilibrium 
on a dynamical timescale. Quasi-equilibrium models of 
compact binaries have been constructed both for neutron 
stars (see, e.g., [|j for corotating and |^ for irrotational 
binary neutron stars) and for black holes 



1 



Even though individual QE models only represent 
"snap-shots" of binaries at a certain separation, it is 
possible to construct the complete adiabatic inspiral to- 
gether with the emitted gravitational wave signal using 
the following scheme For each separation, a QE 

model can be inserted as a matter source into Einstein's 
equations (this is the "hydro-without-hydro" approach 
demonstrated in Q). Numerically integrating Einstein's 
equations for this given matter source will then yield the 
gravitational wave signal and luminosity of the binary at 
that separation. Interpolation between a discrete sample 
of separations then yields the gravitational wave lumi- 
nosity as a continuous function of separation. Combin- 
ing this with the binary's binding energy as a function of 
separation, one can construct the inspiral rate at all sep- 
arations, and hence the separation as a function of time. 
Given a suitable parameterization, the entire continuous 
inspiral wavetrain can then be constructed. The viabil- 
ity of this approach has recently been demonstrated in ||] 
where a prototype calculation was presented for corotat- 
ing binary neutron stars obeying a polytropic equation of 
state. Note that unlike post-Newtonian approaches, the 
QE scheme uses the fully nonlinear Einstein equations. 
The QE approach is also computationally very efhcient, 
since it is possible to produce a very large number of cy- 
cles from a small number of QE configurations, each of 
which needs to be followed for only a couple of orbital 
periods. 

In this paper, we evaluate the QE approach for a model 
problem in relativistic scalar gravitation, and show that 
it produces excellent agreement with an exact numerical 
solution. While scalar gravity (see, e.g., problem 7.1 in 
Misncr, Thorne and Wheeler ||^) is not the correct the- 
ory of gravitation, it is conceptionally much simpler than 
general relativity (GR), and shares many of its character- 
istic features. This makes scalar gravity a very attractive 
framework for calibrating the QE approximation scheme, 
since, unlike in GR, we can directly compare its results 
with the readily producable exact solution. 

In scalar gravity, the gravitational field is described by 
one scalar function. Also, scalar gravity generates grav- 
itational radiation even in spherical symmetry, so that 
the generation and emission of gravitational waves can 
be studied with a spherically symmetric numerical code, 
involving only one spatial coordinate. This enables us to 
track the effect of radiation reaction exactly over many 
dynamical timescales. Moreover, outgoing wave bound- 
ary conditions can be imposed correctly at arbitrarily 
close separations from the sources in spherical symmetry, 
which eliminates the need for large computational grids. 
The theory also admits a local law of energy conservation, 
while GR only obeys global energy conservation. In nu- 
merical work, such a conservation law provides a strong 
check on the accuracy of integration. For all of these rea- 
sons, scalar gravitation has been employed successfully 
to develop many tools of numerical relativity (see, e.g., 
Shapiro and Teukolsky hereafter ST, and also |11|) 
and we extend that tradition here. 



In this paper we study, in the framework of scalar grav- 
itation, the damped oscillations of a spherical matter 
shell, which we adopt as a simple spherically symmet- 
ric analogue to binary inspiral in GR. A relativistic bi- 
nary emits gravitational waves, which slowly extracts en- 
ergy from the binary orbit, so that the inspiral proceeds 
along a sequence of nearly periodic circular orbits of de- 
creasing separation. Similarly, an oscillating matter shell 
emits gravitational waves, which slowly extracts energy 
from the oscillation, so that the damping proceeds along 
a sequence of nearly periodic oscillations of decreasing 
amplitude. Similar to relativistic binaries, where we con- 
sider the quasi-circular orbits at a certain separation of 
the adiabatic inspiral to be in QE, we may also consider 
the matter shell's quasi-periodic oscillations of a certain 
amplitude to be in QE. One difference between the two 
processes is that the binary inspiral continues until co- 
alescence and merger, while the damping of the matter 
shell's oscillations will continuously slow down until a 
true equilibrium state has been reached. 

We adopt three distinct approaches to computing the 
damped oscillation of the matter shell. We first compute 
the exact solution, by numerically integrating the exact 
equations. In our second approach, we adopt a QE ap- 
proximation by neglecting gravitational waves, and com- 
puting QE models of the periodic, undamped oscillations 
of matter shells for various oscillation amplitudes. These 
QE models are then inserted as matter sources into the 
dynamical equations for the gravitational field. Integrat- 
ing the field equations for these QE matter sources yields 
the gravitational wave form and luminosity for each os- 
cillation amplitude. Combining the wave luminosity with 
the QE oscillation energy as a function of amplitude 
yields the amplitude decay rate together with the con- 
tinuous gravitational wavetrain. This is the analogue to 
the QE approach to binary inspiral as outlined above. In 
a third approach, we use a monopole formula to compute 
the gravitational waveforms and wave luminosity for QE 
models and hence to determine the inspiral rate. This 
is the equivalent to using the quadrupole formula in GR 
to compute the inspiral rate and waveforms for binary 
models. 

We compare these three different approaches for three 
different initial configurations, representing weak, moder- 
ately strong and ultra-strong field cases. We find that for 
the weak field configuration, all three approaches agree 
very well. For the moderately strong field configuration, 
the QE approach still agrees very well with the exact so- 
lution, while the agreement with the monopole result is 
much worse. Only for ultra-strong fields, for which the 
assumptions of QE break down, do we find a disagree- 
ment between the QE approach and the exact solution. 
However, even in this case, the break-down is gradual 
and not abrupt. This is a very encouraging result, and 
suggests that the QE approach is a very reliable and effi- 
cient framework for computing adiabatic binary inspiral 
up to the ISCO. 

The paper is organized as follows. We summarize the 
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basic equations of scalar gravity in Sec. II, and outline the 
QE scheme in the Sec. [II. In Sec. [V we present our nu- 
merical results, and compare our three approaches for the 
three different situations. We summarize and discuss the 
implications of our findings in Sec. 0. We also include 
three Appendices. Appendix A contains a short proof 
that a minimum in the quasi-cquilibrium energy corre- 
sponds to a static equilibrium shell solution. Appendix B 
presents some details of the QE scheme for constructing a 
continuous wavetrain. Appendix C describes our numer- 
ical implementation of the field equation. Throughout 
the paper we adopt gravitational units with G = c = 1 . 



II. THE BASIC EQUATIONS 

A. Dynamical equations 

We follow exercise (7.1) in Misner, Thorne and Whee- 
ler and study test particles in a relativistic scalar grav- 
itational theory (see also ST). The field equation for the 
scalar gravitational field $ is 

□ $ = 47re*/9, (2.1) 

where the metric is the flat Minkowski metric 

5o/3 = rfj. (2.2) 



Note that the exponential term on the right hand side 
makes the field equation nonlinear. For a single particle 
of rest mass m, traveling along its worldline z^{t), the 
comoving density p can be written 



1^-9 



Po 
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(2.3) 



where 7 



is the Lorentz factor and po is the density 



in the stationary frame. The particle follows a geodesic 



dr 



+ + = 0, 



(2.4) 



where D denotes covariant differentiation and i 
dz^ /dr is the 4- velocity. 

Conservation of energy-momentum follows from 



V • T = 0, 



(2.5) 



where the stress-energy tensor T'*" can be decomposed 
into a gravitational field part and a matter part: 



T. 



field 



article 



where 



field 



rr-iparticle 



(2.6) 

(2.7) 
(2.8) 



Note that the field equation ( |2.lD can be rewritten 

□ $ = _47rTP^'-"<=i°, (2.9) 

where ^^P^^'^icle _ ^/ii^^^partlclc 

Matter conservation is expressed by the condition 

V-./=0, (2.10) 

where the components of the matter current density are 
j" = 7p and J* = ^pv''. Here is the usual three- 
velocity. 

The above equations can be generalized to a swarm of 
particles by letting 



E 



ruA, 



etc. 



(2.11) 



(see equation ( 2.16) and related paragraph in ST). 

Integrating ( ^.10 ) over all space yields a conserved rest 
mass 



Mo 



J J"d^x = J -fpd^x = J pod^x = const. (2.12) 



Combining the integrals of Eq. (2.5) and Eq. ( 2.10| ) gives 
rise to a conserved total energy at any time t inside a 
sphere of arbitrary radius r centered at the origin: 



(2.13) 



Here Ei is the energy of the gravitational field, including 
a dynamical component, 

Ei^ — f r'^dr' f ($2 + (V$)2)dl], (2.14) 
Stt Jo J ' 

E2 is the particle's kinetic and gravitational binding en- 
ergy. 



Eo 



r' dr' / po(e*7- l)df^, 



(2.15) 



and i?3 is the total outgoing flux of particles and radiation 
across r, integrated over all time. 



E. 



-r- I dt' I dn 



-5-$.t$.^-poi'r(e*7- 1) 
47r 



(2.16) 



Since the total energy is conserved, and since E3 vanishes 
initially, i?tot at all times has to be equal to the sum of 
El and £'2 at t = 0, 
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i?tot(0) = J^r''dr'Jdn[-i^', + (V$)2) 
+Po(e*7-l)]L^, 



(2.17) 



As shown in ST, the total conserved mass-energy M ~ 
J T°°d^x is related to Etot according to M = Mo + Etot- 
The radiative energy flux in the wave zone is 
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47r 



and the total rate of energy emission 
dEficid 



dt 



field 



(2.18) 



(2.19) 



Note that these expressions become exact as r — s- cx). 

In the weak-field, slow-motion limit, the radiation field 
can also be expressed as a multipole expansion, which we 
will use to compare with our QE approximation. Since 
the theory involves a scalar field, the lowest-order contri- 
bution to the radiation arises from the monopole term. 
Using the usual Green's function for a wave equation, 
we follow ST and transform the wave Eq. (2T) into the 
integral form 



$(t,r) 



3 / [e*p]rot 



d-'x 





■g2$ - 


j d^x' 





(2.20) 



where have replaced |x— x'| with r — |x| for large separa- 
tions, and where "ret" means evaluate at retarded time 
t' = t — \x — x'\. For a spherically symmetric density 
distribution, the leading-order radiation field of $ gives 
rise to the monopole formula 



$(i,r) 



47r 
r 



dr'r' 



Po(* 



2 ^ 6 ' 



(2.21) 



This equation is the analogue of the "quadrupole for- 
mula" in general relativity (see also the discussion be- 
low eq. (3.8) in ST). Note again that scalar gravity ad- 
mits gravitational radiation even in spherical symmetry, 
in contrast with GR. 



B. Spherical matter shell 

We now consider a thin, spherical shell of coUisionless 
particles, all of the same rest mass ttia- At every point 
on the shell the particles move isotropically in the plane 
perpendicular to the radius. In an oscillating shell, each 
particle moves about the center in a bound orbit. In the 
Newtonian limit, each orbit is a closed ellipse. 



In spherical symmetry, the geodesic equation (2.4) for 
a particle in the shell becomes 



dR Ur 

dUr _ "0 _ 2'S>^,r 
dt ~ U^R^ ^ ' 



U0 



const, 



(2.22) 

(2.23) 
(2.24) 



where 



,-,0 _ . ■2<S> 



^e^'^+ul + ul/R\ 



and where we have defined 



(2.25) 



(2.26) 



Each particle orbits in a plane, conserving its orbital an- 
gular momentum u^. Note that it is sufficient to inte- 
grate the geodesic equations for one particle, which then 
represents the entire swarm. Note also that for a static 
gravitational field, the particle energy vP is constant. 
The particle mass density is 



p = V— <5(r-i?)<5(0-0^)<5(^-0^) 2 



sm f 



(2.27) 



where {OatiPa) are distributed isotropically on a sphere. 
Inserting equation ( 2.27] ) into equation (2.12) yields 



Mo = ^m^, 

A 



(2.28) 



so that smoothing out the particle distribution in the 
angular direction, we may rewrite the density as a purely 
radial function. 



AuR'^-f 



5{r - i?). 



(2.29) 



where R is the radius of the shell. 

In spherical symmetry, the field equation (2T) can be 
written as 



(2.30) 



Regularity at the origin requires the boundary condition 

= 0, at r = 0, (2.31) 

and we impose an outgoing wave boundary condition at 
the outer boundary. 



(r$),t + (r$),, = 0. 



(2.32) 



The delta function on the right hand side of ( 2.30| ) intro- 
duces a discontinuity in the first space derivative of $. 
Integrating the field equation across the shell yields the 
jump condition 



Mo e2* 



(2.33) 



Note that $ itself is continuous across the shell. In equa- 
tion (2.23), the force term <& ^ then has to be replaced 
by 



(2.34) 



This expression can be found by properly averaging <I> j. 
over an extended shell, and then taking the limit as the 
shell thickness goes to zero. 
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C. static solutions 

For static solutions, in which all particles are in circular 
orbits, the field equation ( 2.3C ) reduces to 



-2C 



47re*p. 



(2.35) 



In vacuum, $ is either constant of is proportional to 



Given the boundary conditions (p.3l|) and (2.32) 



we therefore find solutions of the form 

r<Rs, 
r>Rs, 



<f(r) = -C, 
$(r) = , 



(2.36) 



(2C-C2)i/2 



Mo' 



(2.43) 



Given Mq and the angular momentum u^, this equation 
can be solved for C. The result can then be inserted 



into (|2.39D to yield the radius i?s- 

The oscillations which we consider in this paper are 
damped, and ultimately give rise to a static equilibrium 
state. During the damping process, both the rest mass 
Mq and the angular momentum {t^ are conserved. Given 
any initial non-equilibrium con figuration with Mq and 
{t0, we can therefore use (2.43) to determine the final 
equilibrium configuration to which the oscillating shell 
will ultimately settle down. 



where C and Mc are constants, and where Rs is the 
static equilibrium radius. Note that Mc determines the 
motion of distant particles and gives rise to Kepler's laws, 
so that it can be identified with a "Coulomb" mass, as 
discussed in ST. Unlike in GR, Mq does not agree with 
the tota l cons erved mass-energy M = Mq + Etot defined 
in ( |2T^ ) and (|t|). 

Since <& is continuous across the shell, we have 



Mc 
Rs 



and using the jump condition (2.33) yields 



Mo 

Rs uO ■ 



(2.37) 



(2.38) 



For a circular orbits we have — du^/dt ~ 0, so that 
equation ( ^.23| ) yields 



g-2C 



and therefore 



e^^(l + C/2)i/2. 



Inserting this into ( 2.38| ) yields 
Mc Mo e 



-c 



C = 



Rs Rs (l + C/2)i/2- 



(2.39) 



(2.40) 



(2.41) 



Given Mq/ Rs, equation (2.41) can be solved numerically 
for C or, equivalently, Mc/Rs- In the Newtonian limit, 
where Mq/Rs is small, we find 



Mc = Mq 1 - 



5 Mo 
4 i?^ 



(2.42) 



In the limit of large Mq/Rs, C and Mc/Rs scale with 
log(Mo/i?s), implying that unlike in GR, there is no max- 
imum compaction Mc/Rs- 

We can also combine equations ( 2.39| ) and ( 2.41 ) to 
obtain 



D. Newtonian limit 

In the Newtonian limit, i.e. in the limit of weak fields 
and slow velocities, an analytic solution for a periodi- 
cally oscillating shell can be derived (see Sect. VI of ST). 
Consider a static shell of rest mass Mo and radius 
and reduce all velocities instantaneously by a factor ^ 
at t = 0. The individual particles comprising the shell 
then all move in elliptical orbits with the same period, 
eccentricity and semimajor axes satisfying 



R = Rix{t), 



(2.44) 



where x{t) is given by the usual parametric equations for 
an elliptic orbit: 



X = a(l + ecosu), 

t — — (u -|- esmw), 
27r 



(2.45) 



Here the semimajor axis, eccentricity and period are, re- 
spectively 



1 



P = 27r 



(2.46) 



2i?3 X 



Mo{2~e)\ 

The radial and tangential particle velocities are given by 



X 

Vr — —R, 

X 



V4> = ^- 



Mo 
i?? 



1/2 



(2.47) 



We have used this analytic solution extensively to test 
our code in the Newtonian limit. 



Inserting the analytic solution into Eq. (2.21) and dif- 
ferentiating with respect to time yields the wave ampli- 
tude 
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A = rX 



4 M§ r X 

3 Ri 



(wave zone) 



(2.48) 



The total rate of energy emission can now be found from 
Eq. (ra, 



dE _ 16 
lit ~ T"Rf 



(2.49) 



III. QUASI-EQUILIBRIUM SCHEME 

In the QE approximation, we assume that the orbital 
decay time is much longer than the orbital period. On an 
orbital timescale, the effect of the gravitational radiation 
can then be neglected, and it is reasonable to assume that 
each orbit is determined by the "Coulomb" part of the 
gravitational field rather than the radiative part. This 
suggests that we can neg lect the second time derivative 
in the field equation (2.1), so that the field then satisfies 
the elliptic equation 



47re*pQE. 



(3.1) 



Similarly, QE approximations in GR typically lead to 
elliptic equations for the gravitational field components 
(see, e.g., p^,|3|,^). This approximation greatly simplies 
the problem, since the gravitational field no longer has 
any dynamical degrees of freedom. 

Since no radiation is generated in the QE approxima- 
tion, the system is strictly conservative and there is a 
conserved energy -Eqe- This energy can be derived by 
multiplying equation (2.23) with 2i2,., which yields 



dR 
dt' 



(3.2) 



Here $sh — ^{R), and we have used cons ervat ion of angu- 
lar momentum ii^ — const. From Eqs. (2.34) and ( 2.36| ) 
we also have 



2R' 



Using (2.38), we can now rewrite eq. (B.2) as 



I u'^'^ldR 

i?2^ Mo dt ■ 



(3.3) 



(3.4) 



Expressing the left hand side in terms of u'^ (eq. ( 2.25 ) 
we finally find 



QE 



MqU 



1 



-i?$s[j — const. 



(3.5) 



In order to r elate Eqe to the conserved total energy 
Etot (eq. (2.13)), we evaluate the latt er in the QE ap- 
proximation by inserting the solution ( ^.36 ) and by set- 
ting the radiative components of the field equations to 
zero: $ ff = $ f = 0. We then find 



E^ = ll (v<f)W, 



and 

E2= I r'^dr' I po{u" - l)dn = Mo{u" ~ 1). (3.7) 



The radiative contribution E^ in (2.13) vanishes identi- 
cally. We therefore have 

i?tot =E,+E2= Mo{u^ - 1) - l-'i>l,R^{- - ^), (3.8) 

2 r K 

and find for r ^ cx) 



Etnt — El 



QE 



Mq — const. 



(3.9) 



Note that this energy is conserved only when Ei is evalu- 
ated with r — > oo, so that it includes the entire potential 
energy of the longitudinal (or Coulomb-like) gravitational 
fields. 

In the Newtonian limit. 
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l + $s 



1 



Since to lowest order ~ ^Mq/R, we have 

1 , M2 
= 2^^°^ - 2i 



(3.10) 



(3.11) 



in the Newtonian limit. 

In the bottom panel of Fig. ^ we show Eqe as a func- 
tion of apocenter radius -Rap for constant angular mo- 
mentum U0 = 1.1 7Mq (solid line). As we show in Ap- 
pendix the turning point of this curve corresponds to 
an equilibrium configuration in which all particles move 
in circular orbits (see also Sec. |II C[ ). As the damped 
oscillations of the matter shell radiates energy, the shell 
"slides down" the energy curve in Fig. |l], until it settles 
down to static equilibrium at the curve's minimum. We 
also include in the bottom panel of Fig. |l| the result of 
an exact evolution (dashed line) , in which we started the 
oscillation with a static shell model at Ri — SMq, and 
then suddenly reduced all particle velocities by a factor 
of ^ = 0.7, so that the particle's angular momentum is 
again = 1.17Mq. 

We now adopt a perturbative approach, in which we 
insert the predetermined QE matter density pqe as a 
source to the fully dynamical field equations 



□ $(r, i) = 47re*pQE, 



(3.12) 



where pqe is given by eqn. (2.29) using the QE solu- 
tion for R{t) (compare with the "hydro-without-hydro" 
approach as suggested in Q). Integrating this equation 
then yields the periodic gravitational wave form and lu- 
minosity dE/dt for each QE configuration of a certain 



6 



<dE/dt> 




0.947 



FIG. 1. The gravitational wave luminosity (upper panel) 
and the total energy (lower panel) as functions of apocenter 
radius ( for a m oderately strong field configuration, Ri = 8Mo, 
below). The dotted line denotes the result from 
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see Sec. 

the monopole formula (2.49), the dashed line shows the QE 



result, and the solid line shows the integration of the ex- 
act equations. I n th e lower panel, the dashed line is com- 
puted from eq. (3.5), and the solid line is computed from 
-Eqe El + E2 + Mo as in eq. (^.13|) with r ^ 00. 



apocenter radius. In practice, we determine dE/dt by 
averaging over one period P 



dE 



1 



< ^ >= - / ^dt' 



dt 



P 



dE 



dt 



(3.13) 



We show < dE/dt > as computed from QE models in the 
top panel of Fig. |l] (soUd line) and compare both with re- 
sults from the integration of the exact equations (dashed 
line) and the monopole formulae (dotted line) . From this 
plot it is already obvious that the QE "hydro-without- 
hydro" approach provides a much better approximation 
than the monopole formula. Combining < dE/dt > for 
several values of i?ap with the derivative of the QE energy 
dEQE / dRiip then yields the damping rate 



dR, 



dt 



< dE/dt > 
dEQE./dRa.p 



(3.14) 



(This is analogous to eq. (1) in Q for binary neutron 
star inspiral.) Note that at the final equilibrium ra- 



dius R = Rs, the numerator < dE/dt > vanishes 
more rapidly then the denominator dEq^ / dR^p , so that 
dRa^/dt smoothly approaches zero. 

The complete wavetrain of the damped oscillations can 
be assembled by suitably parameterizing the wave sig- 
nals for different values of i?ap (compare with eq. (2) 
in [^). We p rovide details of this parameterization in 
Appendix B 2 . 



IV. NUMERICAL RESULTS 

In this section, we compare results from the QE ap- 
proach, the exact integration and the monopole formula 
for weak, moderately strong and ultra-strong fields. For 



each case, we prepare a static shell solution as in Sec. [I C 



for a certain value of Mq / Ri , where Ri is the initial apoc- 
enter radius, and then reduce all velocities and hence an- 
gular momenta by a factor of ^ = 0.7. The particle orbits 
are then out of equilibrium and start a damped oscil- 
lation, until they lose sufficient energy by gravitational 
radiation to settle down into the final static equilibrium 
corresponding to the reduced value of the angular mo- 
mentum. 

Given the above scenario, the appropriate initial con- 
ditions for the gravitational field are 



0, 



(4.1) 



where $static is given by Eqs. ( p^ ) and ( ^Isj) . The 
particle initial data are 



R — Rs, 



iir = 0, 



and ua, = £,uZ 



(4.2) 



where {i^tatic jg gjygjj j^y ( |2.37| ) and (|2.39|) . Thus the (non- 
equilibrium) shell begins at rest with all the particles at 
their apocenter positions. 

For numerical reasons (involving the regridding algo- 
rithm as described in Appendix y), we found it conve- 
nient to impose outer boundary conditions at 5Ri (except 
for the strong field case, where we choose 50Ri for the 
outer boundary). We resolve the interior region of the 
shell with 50 gridpoints, and the exterior region with 200 
gridpoints. The field integration scheme is adopted from 
ST and summarized in Appendix 



A. Weak field configuration {Ri = lOOOMo) 

We first study a weak field configuration with Ri = 
IOOOA/q. After having reduced the angular momentum 
be a factor of ^ = 0.7, this configuration will ultimately 
settle down to a final radius of Rf ^ Rs ~ 491. BMq. 
In Fig. H we compare the evolution of the shell's ra- 
dius as fou nd from the integration of the exact equations 
(Sec. |II B[ ), the QE approach (Sec. Ill), and the analytic 
Newtonian result (Sec. II D| ). As expected, the agreement 
between the different approaches is excellent. 
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FIG. 2. Evolution of the radius as a function of time for a 
matter shell with the initial radius Ri = lOOOMo and velocity 
cut-down factor ^ = 0.7. The solid line is the result from the 
integration of the exact equations. The dotted line, which 
completely coincides with the solid line, marks the analytic 
Newtonian result of Sec. [ID. The short-dashed line is the 



envelope of apocenter radii oscillations according to the QE 
approach. The long-dashed line marks the final equilibrium 
radius of the shell. 



FIG. 3. Contributions to the total energy as a function 
of time for the weak field case Ri = lOOO Mp. The solid line 
shows the total energy i?tot (see eq. (2.13). The field energy 



El is marked by the dotted line, the particle energy E2 by the 
long-dashed line, and the integrated flux by the dot-dashed 
line (see also top panel). We also include the initial value of 
-Etot(O), marked by a horizontal dashed line. The agreement 
between the solid and the dashed line demonstrates energy 
conservation and is a measure of the accuracy of our code. 



In Fig. 1^ we show that energy is conserved to very high 
accuracy in o ur code. Here, we evaluate the conserved 
integral ( 2.17 ) at radius r = IGOOAfg. The plot also shows 
that the radiated energy E^, (integrated flux) is very small 
compared with the other energy terms. 

We show the QE parameters from which the QE wave- 
form is constructed in Fig. ^ The gravitational radia- 
tion luminosity < dEq^/dt > is computed by integrating 
eq. ( 3.12| ) for 13 different values of j?a p and interpolat- 
ing between the results. Using ( |3.13 ) and ( |3.14 ), this 
can be combined with dEq^/ dR^^ to yield the damping 
The 4 parameters ea, P, and ep of the 



rate dRi,p/dt_ 



are 
As 



waveform (B12) and (B13) at different values of i?ap 
determined by nonlinear data fitting (Appendix ^). 
expected, the values of ca and ep are very similar in the 
weak field case. 

In Fig. H we compare the waveforms as obtained from 
the integration of the exact equations (solid line), the 
QE approach (dashed line) and the Newtonian analytic 
solution (dotted line) for the first 30 periods. The QE 
approach can reproduce the exact result very well. Note 
also that the QE scheme is very efficient: given an in- 
terpolation between the data for a small set of apocenter 
radii i?ap, the entire wavetrain, which in this case would 
take thousands of cycles and would, in an exact integra- 
tion, be dominated by accumulation of numerical noise, 
can be constructed very easily. 



B. Moderately strong field configuration {Ri — 8A/o) 

We now turn to a configuration with a moderately 
strong field, Ri = SA/q. In this case, the oscillation 
damps out much more quickly due to radiation reac- 
tion, and we can integrate the exact equations until equi- 
librium has essentially been reached. This occurs at 
R = 4.77Mo. 

We show the evolution of the shell's radius in Fig. ||. 
This plot includes the result from the integration of the 
exact equations (solid line) as well as the "envelope" 
Rap{t) as found in the QE approach. Fig. shows the QE 
parameters which we have constructed for this configu- 
ration. We plot the energy contributions as a function 
of time in Fig. ^, and find that energy is conserved to 
about 0.4%. Computing the exact solution takes only a 
few CPU hours on, for example, an SGI 02 workstation. 

Fig. ^ shows three waveforms of the oscillating shell. 
Here we compare the results of the QE and exact inte- 
grations with the result of the monopole formula when 
applied to the oscillating QE configurations. The QE 
waveform and the exact waveform agree very well up to 
late times {t ~ 7000Mo), at which point the oscillation 
amplitude has decreased by about a factor of 100. At 
this point the integration of the exact equations has ac- 
cumulated substantial numerical noise, and may actually 
be less accurate than the QE approach. The waveform 
from the monopole approach disagrees with the exact one 
even at very early times, showing that the QE approach 
is much more reliable for moderately strong fields. 
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FIG. 4. Calculated QE data and their fitting curves as 
functions of apocenter radius for the Ri — lOOOAfo case. The 
circles are the calculated data points, and the dotted lines are 
the fitted curves. 



C. Ultra-Strong field configuration (R, = O.lMo) 

For very strong fields, the oscillations are damped so 
strongly that our QE assumption no longer holds, and 
we therefore expect our QE approach to break down (see 
discussion below). Fig. shows that the QE envelope 
Ra.p{t) no longer matches the exact result very well. We 
show the fitting parameters for this case in Fig. 0, and 
the energy contributions in Fig. 12 which still obey en- 
ergy conservation quite well. 

In Fig. |l^ we compare the waveforms from the ex- 
act integration and the QE approach, and find the ex- 
pected disagreement. Even here, however, the overall 
shape of the QE waveform is not far off. We do not in- 
clude the monopole approximation, since its predictions 
do not even fit on the same scale. 

We can illustrate the breakdown of the QE approach 
by the following simple argument. Figures |, I and |ll| 



show that dRiip/dt decays like 



dt 



-a(i?ap - Rf), 



(4.3) 



where the value of a depends on Ri/Mo and can be read 
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FIG. 5. Wave amplitude A (multiplied by r) as a function 
of time at r = IGOOAfo for the weak field case {Ri = lOOOAfo). 
The solid line shows the integration of the exact equations, 
the dashed line marks the QE result, and the dotted line is 
the monopole radi ation obtained from the analytic Newtonian 
solution (see Sec. II D). 



off of the above figures. As a result, i?ap decreases expo- 
nentially. 



f 



Rj — R 



f 



(4.4) 



The timescale trad = ce^^ can now be compared to the 
orbital timescale iorb- Inserting numbers for the three 
cases, we find 



^rad 
^oib 



280000, 


Ri 


= lOOOA/o 


15, 


Ri 


= 8A^o, 


0.5, 


Ri 


= O.lAfo. 



(4.5) 



For the QE assumptions to hold we need Q ^ 1, which 
obviously holds only for Ri = lOOOAfo and iJ^ = SAfo, 
but not for the ultra-strong field case Ri — 0.1 Afo- It 
is therefore not surprising that the QE approach breaks 
down in this case. 

We may also check the QE a ppr oximation, which ne- 
glects $.ft in the field equation (2T), for self-consistency. 
Assuming that $ ~ 1/i?, we estimate its magnitude to 
be 



R 



2R^ 



R^ 
i?3 



and similarly 



- — 



1 



(4.6) 



(4.7) 
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FIG. 6. Evolution of the shell radius as a function of time 
for a matter shell for the moderate field case {Ri = 8Mo). 
Labeling is the same as in Fig. ^ except that we do not include 
the Newtonian result. 



The ratio of the two terms in equation (|2.lD then becomes 



V2$ 



vl. 



(4.8) 



This suggests that $,ft can be neglected, and hence that 
the QE approximation can be applied, only when <C 
1. In Figure |lj we show the maximum radial velocity 
for different values of Ri/Mg, and also plot the ratio of 
^ tt and 2^ r/R (which scales with V^<I>, but is finite at 
the shell). This plot suggests that for Ri/Mo > 1 we 
should expect errors less than a few percent, while for 
Ri/Mo = 1, as in our ultra-strong field case, we should 
expect errors on the order of 10 %, which is consistent 
with our results. 



V. CONCLUSION 

In this paper, we illustrate our QE approach for cal- 
culating gravitational waveforms for a model problem in 
scalar gravitation. We demonstrate why it is a very vi- 
able technique to modeling binary inspiral in GR. 

Compact binaries, consisting of neutron stars or black 
holes, emit gravitational wave and slowly spiral towards 
each other until they reach the ISCO. Outside of the 
ISCO, the inspiral is very slow: the gravitational radi- 
ation reaction timescale is much longer than the orbital 
period. It is therefore reasonable to assume the bina- 
ries to be in QE, and the inspiral to proceed along a 
sequence of QE configurations. In a recent paper, Duez, 
Baumgarte and Shapiro |^ demonstrated how the inspi- 
ral can be modeled by inserting QE binary models as 
matter sources in Einstein's field equations. Integrating 
Einstein's equations then yields the gravitational wave 
luminosity, from which the inspiral rate and hence the 
entire gravitational wavetrain can be constructed. 
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FIG. 7. Calculated data QE and their fitting curves as 
functions of apocenter radius for the moderately strong field 
case {Ri = 8Afo). 



In this paper, we study a simple analogue in scalar 
gravitation. Scalar gravitation has the advantage that it 
is conceptionally simpler, and that it admits gravitational 
radiation even in spherical symmetry, so that we can re- 
produce all qualitative characteristics of the QE approx- 
imation even in a 1 + 1 dimensional problem. We study 
the radiative damping of an oscillating, spherical matter 
shell. In analogy with the orbiting binary in GR, the 
oscillating shell emits gravitational radiation, and hence 
looses energy. While the binary separation decreases in 
the process, the oscillation amplitude decreases. Ulti- 
mately, the oscillation is completely damped out, and a 
static equilibrium configuration is reached. 

The advantage of scalar gravitation is that it is pos- 
sible to integrate the exact equations without any ap- 
proximation. We compare these exact results with pre- 
dictions from both our QE approach and the monopole 
approximation (which is the analogue of the quadrupole 
approximation in GR) . We find that all three approaches 
agree very well for weak field solutions, but that the QE 
approach reproduces the exact solution much better for 
moderately strong fields. Only for ultra-strong fields does 
the QE approach break down. 

We conclude that the QE approximation is a very 
promising approach to computing the adiabatic inspiral 
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FIG. 8. Energy conservation as a function of time for the 
moderately strong field case {Ri — 8Afo). Labeling is the 
same as in Fig. pi 
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APPENDIX A: THE CORRESPONDENCE 
BETWEEN A TURNING POINT IN THE 
ENERGY CURVE AND STATIC EQUILIBRIUM 

In this section we show that a turning point of the 
quasi-equilibrium energy E'qe versus apocenter radius 



QE 



MqU 



1 



(Al) 



along a sequence of constant rest mass and angular mo- 
mentum, coincides with an equilibrium configuration in 
which all particles are in circular orbits of radius Rc- 

We e valua te -Eqe at apocenter radius R^p, where Ur = 
0. Eq. (|2.25| ) then yields 



(A2) 



For a sequence of constant rest mass and angular mo- 
mentum we therefore find 



d(Mou°) 

dRan 



-i?ap$sh 



dRav 



uHloRlp ' 



(A3) 



FIG. 9. Wave amplitude A (multiplied by r) as a function 
of time as measured at r = lOMo for the moderate field case 
{Ri — SMo). Labeling is the same as in Fig. ^, except that 
the dotted line denotes the monopole formula applied to the 
oscillating QE shell. 



where we have substituted eq. (2.38) in the form 



u°R 



(A4) 



ap 



to get eqn. (A3). The derivative of the second term in 
eq. (P) is 



dR^^^ - ^^^^'""dK: 



Combining eqs. (A3) and (A5) we now find 



dE^_^_Mo^ 
dR: 



ap 



(A5) 



(A6) 



At a turning point dEQ^jdR^-p — 0, which implies 



-Rl^sh- 



(A7) 



Equation ( A7) is equ ivalent to the result for static spher- 
ical shells, eq. ( ^.39| ). A turning point in the quasi-equi- 
librium energy therefore implies 



^ap — 7 



(A8) 



and hence coincides with a static equilibrium configura- 
tion. 
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FIG. 10. Evolution of the radius as a function of time for 
the ultra-strong field case {Ri = O.lMo). Labeling is the same 
as in Fig. ^, except that we do not include the Newtonian 
result. 



APPENDIX B: CONSTRUCTING THE 
CONTINUOUS QE WAVETRAIN 

In this appendix we discuss the assembly of the QE 
solution and waveform. To do so, we construct QE so- 
lutions for a set of apocenter radii. We then evolve the 
gravitational fields dyn amica lly in the presence of a QE 
matter source (see eq. ( ^.12[ )) to determine the gravita- 
tional wave form and luminosity for each apocenter ra- 
dius. Guided by the analytical Newto nian expression for 
the shape of the waveform (see Sec. 11 D), we parame- 
terize the waveform to be able to match the computed 
periodic data at discrete radii to a smooth function. In 
this Appendix we describe this parameterization together 
with a prescription for how the time evolution of these pa- 
rameter s ca n be determined. We illustrate our approach 
in Sec. 



Bl 



where we focus on the phase of the gravi- 
tational waves and its dependence on other par amet ers, 
and we construct the waveform template in Sec. B2. 



1. Waveform phase versus time 

Consider the simple example of an oscillator with con- 
stant period P and constant wave amplitude A. The 
waveform A can then be found from the phase 9, which 
in turn is given as a function of time: 



A 



2n 

Asin( 



(Bl) 
(B2) 



In general, however, both A and P vary with time, in 
which case the phase, now denoted by 'd can be computed 
from an integration 
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FIG. 11. Calculated QE data and their fitting curves as 
functions of apocenter radius for the ultra-strong field case 
{Ri = O.lMo). Labeling is the same as in Fig. ^ 
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(B3) 



In some cases, it is convenient to introduce a phase 
angle such that 

d = d + ip, (B4) 



where 9 is given by eq. (Bl). Taking a time derivative of 
this equation, we find 

27r , ■ 2tt 2Trt ■ 



P 



P P2 



(B5) 



where the first equality follows from (B3), the second 
from (B4), and the third from (Bl). Comparing the left 
and right hand side, we find (p = 2'KPt/P^. 

The relation between the phase 9 and time t may 
be more complicated than eq. (Bl), as for example in 



the Newt onian analytic solution for an oscillating shell, 
eq. ( 2.45 ). We therefore allow = 9{t,pi) to be a func- 
tion of time t and a set of additional parameters Pi(t), 
and generalize eqs. ( p3|) to 



J dt9dt, 



(B6) 
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FIG. 12. Energy conservation as a function of time for the 
strong field case {Ri = O.lMo). Labeling is the same as in 
Fig. I 



We again introduce the phase angle (p 

d = e + ^, (B7) 
dtB = d = e + Lp = dte + i:p,dp£ + Lp, (bs) 



and now find 



(B9) 



This expression determines how the phase angle evolves 
when 9 depends on several parameters. We will use this 
expression for the scalar wave model problem below. 



QE parameterization and construction of the 
wave template 



As described in Sec. Ill, we construct QE models for a 
set of apocenter radii and insert these as matter sources 
into the field equations as matter sources (eq. ( 3.12 )). 
Dynamically evolving the gravitational fields will then 
yield the gravitational wave luminosity and the gravita- 
tional wave form. Given a suitable parameterization of 
the wave form, whose shape may be far from sinusoidal, 
we can find a set of parameters for each chosen value of 
the apocenter, and can then construct an interpolating 
function which yields all parameters as smooth functions 
of the apocenter (e.g. Figs. H, |^ and |ll|). Combining 
the QE energy with the gravitational wave luminosity 
(eq. (3.14)) yields the damping rate and the entire evolu- 
tion of the system. This enables us to express all param- 
eters, and hence the waveform, as a continuous function 
of time. 

A suitable wave from template can be con struc ted from 
the New tonia n analytical solution (see Sec. [I D ). Insert- 
ing eq. (|2.45|) into (|2.48|) yields 
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FIG. 13. Wave amplitude A (multiplied by r) as a function 
of time as measured at r = O.8M0 for the ultra-strong field 
case [Ri = O.lMo). Labeling is the same as in Fig. ^, except 
that we do not include the Newtonian result. The small bump 
at early times is a numerical artifact and is due to slightly 
imprecise initial data. 



t= —{u + e sinu), 
27r 



A: 



Stt M, 



2 r 



ap 



esmu 



aP{l + ecosw)^ 



(BIO) 
(Bll) 



We found that the fitting of the computed waveform to 
the above function in strong field configurations could 
be improved by allow ing the eccentricity e in the two 
equations in eq. ( 2.45 ) to be different. This yields the 
template 



t — — (u + ep sm u) 
2tt 



A = 



A sinu 



(1 -f cosu)'^ 



(B12) 
(B13) 



We thus have four independent parameters P, ep, A, ca 
which we choose to characterize the waveform. 

We construct QE configurations for 15 different apoc- 
enter radii i?ap between the initial radius Ri and the final 
circular radius Rf. At each value of i?ap we determine 
the four parameters P, ep, A, ca as well as the gravi- 
tational wave luminosity, AE/P. We then construct in- 
terpolating functions, so that these parameters are now 
given as smooth functions of i?ap- Given the wave lu- 
minosity AE/P and the QE energy Eq e, th e damping 
rate dRap/dt can be computed from eq. ( [3.14| ). The four 
parameters P, ep. A, ca, the wave luminosity AE/P 
and the damping rate dRg^p/dt are now all available as a 
function of apocenter radius -Rap- We show our results 
for weak, moderately strong and ultra-strong field cases 
in Figs. |, I and ^ 

As in Sec. |B l|, we allow for a phase shift by introducing 
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where u satisfies eq. ( p312| ), and where, according to 
eq. (B9), the phase angle (f evolves according to 



if : 



-epdepU — Pdpu 



2TrtP + P'^ep sinu 



P^(l + ep cos u) 
The wave form A is now given in terms of by 
Asm'3 



A: 



(1 + eA cos-dy 



(B15) 



(B16) 



Here, the time derivative of the four parameters can be 
derived from the chain rule 



Pi 



dpi dR. 
dRap dt 



ap 



(B17) 



The continuous waveform A can now be constructed by 
integrating the two equations 



ip = 



2tt — uP — {epP + Pep) sin u 

P{l + epcosu) 
2TrtP + P'^ep sinu 
P2(l + ep cos u) 



(B18) 

(B19) 
(B20) 



The sum of the two yields whic h, to gether with A and 
eA, can then be inserted into eq. (B16) to give the entire 
gravitational wave train. 



APPENDIX C: NUMERICAL ALGORITHM FOR 
INTEGRATING THE SCALAR WAVE 
EQUATION 

In this Appendix, we outline our finite difference 
scheme, which is based on that of Ref. [|lO|. and ex- 



plain how the jump condition ( 2.33 ) is implemented in 
our code. 



In spherical symmetry, the field Eq. (|2.9| ) is 
1 



(CI) 



We rewrite this second order equation as two equations 
which are first order in time 



FIG. 14. The maximal value of the square of radial velocity where 
Vr (circles linked by the dotted line) and the dynamical part 
of the field equation J?"I>,tt/2$,r (squares linked by the dashed 
line) vs the initial radius. The size of these terms is a rough 
measure of the expected relative error in the QE approach. 



r[$] = A, 

T[A] = 7^[$] + 4ttT, 



^ Yt, 

n[Y] = epr^s]^^ 



(C2) 



(C3) 
(C4) 



Here Y denotes either $ or A. The Laplacian in the 
field equation is written in the form of ( |C4| ) to ensure 
regularity of the finite-difference operator near the origin. 
We implement a leapfrog scheme with a variable time 
step to solve these equations, and finite difference the 
operators T[Y] and TZ[Y] according to 
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where Atn = tn+i —t„ . These operators are second-order 
accurate in both space and time. At r = Tmax we impose 
an outgoing wave boundary condition 



irY),t + (ry),, = 0, 



(C7) 



where Y is either $ or A. A second-order accurate finite 
difference form of this equation is 



Y. 
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(C8) 



where 
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c = 



,c + l/2 



-1/2 



(C9) 



The jump condition (2.33) is easiest implemented by 
letting the grid move with the matter shell , so t hat i?sh — 
ri^j^_i_i/2 at all times. The jump condition ( |2.33| ) can then 
be discretized to 
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i3h + l/2 



i3h-l/2 



^ish+3/2 - ?'i,h + l/2 
A^O 2*" ,,,, 



''j.h + 1/2 



''i,h-l/2 

(CIO) 
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Eq. (CIO) is used (by using a rootfinder) to derive the 
value of $"■"'"1, ,„ with the known values of and 

JBh + l/2 iBh+3 /2 

'^tiich are obtained from eqs. ( |C5| ) and ( |C^ ) in 
advance. After each time step, we interpolate the func- 
tion values (at all three timelevels n + 1, n and n — 1) to 
a new grid, so that 7'igjj+i/2 at the new time level n + 1 
coincides with the new location of the matter shell i?sh- 
In our runs, we have chosen igh = 50. We use a uniform 
grid inside the shell, and a geometric progression, i.e., 
ri+3/2 - ri+i/2 = a{ri+i/2 - ri_i/2) where a is a fixed 
ratio, near unity, in the exterior. 

Since all the particles making up the shell are identical, 
it is sufficient to evolve one particle. M oreov er, sin ce 
is a constant of the motion, only eqs. (2.22) and ( 2.23| ) 
need to be integrated to determine the particle's geodesic 
motion. 
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